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Conventional models of Josephson junction dynamics rely on the absence of low energy quasipar- 
ticle states due to a large superconducting gap. With this assumption the quasiparticle degrees of 
freedom become "frozen out" and the phase difference becomes the only free variable, acting as a 
fictitious particle in a local in time Josephson potential related to the adiabatic and non-dissipative 
supercurrent across the junction. In this article we develop a general framework to incorporate the 
effects of low energy quasiparticles interacting non-adiabatically with the phase degree of freedom. 
Such quasiparticle states exist generically in constriction type junctions with high transparency 
channels or resonant states, as well as in junctions of unconventional superconductors. Further- 
more, recent experiments have revealed the existence of spurious low energy in-gap states in tunnel 
junctions of conventional superconductors - a system for which the adiabatic assumption typically 
is assumed to hold. We show that the resonant interaction with such low energy states rather than 
the Josephson potential defines nonlinear Josephson dynamics at small amplitudes. 

PACS numbers: 74.50+r, 74.78.Na, 72.10.Bg 



I. INTRODUCTION 

During the last twenty years the microscopic theory of 
the Josephson effect has been undergoing steady develop- 
ment following the advent of novel mesoscopic Josephson 
structures such as transparent metallic and semiconduct- 
ing junctions [T], quantum point contacts [2], quantum 
dot contacts [3], junctions with spin-active interfaces [J. 
Much of the theory development for these structures were 
based on pioneering work by I.O. Kulik [5-8 . Also im- 
portant breakthrough was experimental demonstration 
[SHU] of macroscopic quantum coherence [T2] in Joseph- 
son junctions, and realization of quantum Josephson cir- 
cuits (qubits) P3HH>]. 

Functioning of quantum Josephson circuits is based 
on a fundamental property of Josephson tunnel junc- 
tions: nonlinear non-dissipative phase dynamics. Equiv- 
alence of Josephson junctions to ideal nonlinear oscilla- 
tors, pointed out already by Josephson [T7], is used in 
numerous applications in microwave electronics |18j . The 
possibility to quantize the motion of Josephson oscillator 
|19j . and to observe the macroscopic quantum dynamics 
is essentially based on this fundamental property. 

Equation of motion for the superconducting phase dif- 
ference across the junction stems from Kirchhoff 's rule 
that combines the Josephson tunneling current, Ij(tp) = 
Ic sin ip, and the displacement current through junction 
capacitor, (C/2e)<p, 
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where I e (y>, t) is a biasing current defined by external cir- 
cuit, h = 1. A key assumption behind this equation is a 
quasi-static form of the Josephson current that extends 
the static current-phase relation to the non-stationary 
case of temporal variation of the phase. A justification 
for this assumption is provided by a wide isotropic su- 
perconducting energy gap A that prevents excitation of 
quasiparticles by temporal variation of the phase at low 
temperature and small frequency of Josephson plasma 
oscillation, kT, huj p -C A. Thus electrons in the junc- 
tion remain in equilibrium, and the adiabatic form of the 
Josephson current is maintained. 

Such an approach is relevant for tunnel junctions, but 
it is not always correct. Notable exceptions are transpar- 
ent point contacts [2U] and resonant quantum dot con- 
tacts [2 1 J containing Andreev bound states deep inside 
the energy gap. Other important exceptions are junc- 
tions of d-wave superconductors containing zero energy 
Andreev surface states [22] and low energy nodal quasi- 
particles [23j . In such junctions the low energy quasipar- 
ticles are involved in the macroscopic dynamics: they are 
excited and driven away from equilibrium by temporal 
variation of the phase resulting in significant modification 
of the Josephson current. How is Eq. ([I]) then modified 
in the presence of low energy quasiparticle states? 

In this article we suggest an extension of Eq. (JlJ to 
describe the non-adiabatic Josephson dynamics in the 
presence of interaction with quasiparticles. A general 
equation derived in the next sections has the form, 



ip + I c smip = I e (tp,t), 



(1) 
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<p + Tr(Ijf)=I e , 
i 'f = [H, /]• 



(2) 
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Here the adiabatic Josephson current is replaced by a 
statistical average of a Josephson current operator, Ij; 



2 



the non-equilibrium quasiparticle density matrix / satis- 
fies the Liouville equation with an effective Hamiltonian, 
H. The only approximation made during the derivation 
is a semiclassical approximation for the phase dynamics, 
otherwise this is an exact equation. As we will show, 
both the current operator and effective Hamiltonian are 
expressed through the quasiparticle energy spectrum of 
the junction and interlevel transition matrix elements. 

Eq. ([2]) has a generic form of equation of motion of a 
macroscopic particle interacting with a fermionic bath. 
Usually such problems are treated assuming an equilib- 
rium bath. Here we will consider a non-equilibrium bath 
consisting of low energy bound Andreev states strongly 
driven by the phase dynamics. Our main conclusion is 
that the Rabi dynamics of the Andreev states dramat- 
ically modifies the nonlinear properties of macroscopic 
Josephson dynamics. The physics here resembles well 
known in nonlinear optics picture of interaction of elec- 
tromagnetic mode with medium of two- level atoms [24j . 

The structure of the paper is as follows. In section II 
we discuss a general approach based on the path integral 
technique, which is used in Section III to derive Eq. ([2]). 
In the next section we discuss the adiabatic limit and 
establish connection between our method and earlier re- 
sults for tunnel junctions. Section V is devoted to non- 
adiabatic effects; we study both the linear and nonlinear 
quasiparticle response, the main result here is the evalu- 
ation of a nonlinear effect of driven low energy Andreev 
bound states. In Section VI we present the derivation of 
stochastic Langevin equation generalizing the determin- 
istic Eq. 



II. FORMULATION OF THE PROBLEM 

Consider a general setup of a junction with super- 
conducting electrodes occupying left (x < 0) and right 
(x > 0) halfspaces, with an interface at x = carrying N 
conducting modes. We will not specify the properties of 
the interface but rather characterize it, within the qua- 
siclassical approximation, with some electronic transfer 
matrix. In the following we also adopt common assump- 
tions: (i) the superconductors are described with a BCS 
mean field theory, (ii) superconducting electrodes main- 
tain local equilibrium implying absence of spatial and 
temporal variation of the module and phase of the order 
parameter, A = const, x(r, t) = sign(x)<p(i)/2. 

To accurately describe the nonequilibrium dynamics 
we adopt the path integral approach, introduced by Am- 
begaokar et al. [55], and adapted for non-equilibrium 
systems J2()-28|. Following this approach we represent 
the trace of the time dependent statistical operator of 
the junction p(t) = U(t,to)poU(to,t) with the path inte- 
gral, 

Z = Tr[p(£)\ = I VtpVipV?Pe lS[v> ^^ ] , (3) 



where the action is, 

SM,<I>] = Jdt - U B (<p) + ■ (4) 

c 

The first term in this equation originates from the elec- 
trostatic interaction between electrodes and is described 
within the capacitance approximation |25) ; the second 
term is an inductive energy of the external circuit, and 
the last term represents the contribution of supercon- 
ducting electrons. Time integration goes along the 
forward-backward time contour, C = C+ + C_. The 
fermionic fields in the electronic term are written in the 
Nambu pseudo-spinor representation, ip — (ifo , -0^) T , and 

Q- X =id t -U(<p)-\G Z , (5) 

where 

U(<p) = (^--E F + V{T))a z + Ae^a x , (6) 
\2m ) 

is the junction Hamiltonian. Here V(r) is the potential 
defining the interface; superconducting order parameter, 
A, is a scalar in s-wave superconductors, but becomes a 
nonlocal operator in the case of unconventional d-wave 
pairing. The last term in Eq. ^ represents the electrical 
potential needed to preserve electro-neutrality within the 
electrodes [2"5] . 

By virtue of the quadratic form of the Fermionic part 
of the action Q one can formally perform the gaussian 
path integral over the fermionic fields, and reduce the 
integral to one over the phase degree of freedom , 

Z = Jvip e^oM + SplnH^ -1 ), (7) 

here So comprises the first two terms in Eq. Q, and Sp 
denotes the trace over both the quasiparticle states as 
well as the forward-backward time contour. This trans- 
formation in itself, however, does not solve the problem: 
the obtained effective action contains the contour ordered 
fermionic Green's function, which needs to be computed 
by solving the equation of motion. This can only be done 
under some approximations. The most studied in liter- 
ature case concerns tunnel junctions where the Green's 
function is calculated perturbatively using small trans- 
parency of the junction, D <§C 1 [251 I3"T]. This is com- 
monly done within the formalism of tunnel Hamiltonian 
model. This method can be improved and made suitable 
for transparent junctions, D ~ \ } by performing sum- 
mation of the whole perturbative series [34] . However, 
the tunnel model method does not straightforwardly ap- 
ply to superconductors with surface states, such as d- 
wave superconductors, since it is based on expansion over 
bulk Green's functions. The tunnel model must then be 
modified by considering semi-infinite leads with hard- wall 
boundaries rather than homogenous leads |35| . An alter- 
native way to calculate the effective action for transpar- 
ent junctions was suggested in Refs. [2H by using 
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exact boundary conditions and an adiabatic approxima- 
tion for low energy Andreev states. Zaikin and Panuykov 
[521 [33] suggested a general method for calculating the ef- 
fective action by establishing a formal relation between 
the action and the current across the junction. This 
method, however, requires knowledge of the ac current 
response to an arbitrary time dependent realization of 
ip(t), which in general is not possible to obtain. 

In this paper we suggest an alternative method of cal- 
culation of the effective action Q , which is exact in the 
limit of semiclassical phase dynamics, and universal re- 
garding interaction with any kind of quasiparticle states. 



A. Instantaneous Basis 

The central idea of the method is to expand the Nambu 
fields over an instantaneous eigenbasis of the Hamiltonian 



(8) 



This allows us to separate the spatial problem from 
the temporal one by solving the time independent 
Bogoliubov-de Gennes equation for a fixed value of the 
phase. Apart from the technical simplifications this ba- 
sis provides an intuitive understanding of the microscopic 
processes involved in the Josephson dynamics in terms of 
transitions between quasiparticle states. 
In this basis the action Q becomes, 

S , [y,{a i },{a i }]=y alt ^^<P 2 - Uei^+Y^OiG^aj j , 



where 



G- 1 = (id t - H K 



(9) 



(10) 



represents the quasiparticle Green function in the instan- 
taneous basis, and the Hamiltonian is given by equation, 



Hij(<p, <p) = Ei(tp)6ij - ipAi 



(11) 



The diagonal elements here are given by the instanta- 
neous eigen energies of the Hamiltonian ([6|, 



U<\> t = Ei 



(12) 



and the off-diagonal elements are proportional to the ma- 
trix elements, 

Aij = (fa, id v (f>j) - I (fc, sign(x)a z (j) j ) , (13) 

of the transitions between the instantaneous eigenstates 
due to temporal variations of the phase. 

The physical meaning of the transition matrix elements 
can be understood by establishing their connection to the 
Josephson current operator. Consider a general quantum 



mechanical equation for the charge current density ma- 
trix, 



(V-V')^l(r)^(r') 



2m 



(14) 



The current through the interface, S, is given by equa- 
tion, 



dn-jtf(r). 



This is the matrix of the Josephson current operator. If 
we connect the electrodes in a loop at infinity, we can use 
the fact that no current is flowing through any other part 
of the surface of the superconductor so we may extend 
the surface, S, around the whole superconductor and use 
Gauss law: 



2I« = J d 3 r V • jy(r) - J d 3 r V • j tf (r). 



(15) 



From the explicit form of the Hamiltonian ^ we derive 
the relation, 



Hv) h 



e 
2m 



(16) 



The last term in this equation can be rewritten as 

[a z ,H] = 4isign(x) d v H. (17) 
The current operator then becomes 



2ie 



(<j)i,id tp 'U<j)j - ^(Ej - E i )(cj) i ,sign(x)a z (j)j) 



(18) 



By differentiating the eigenvalue equation, H^i = Eiq 
with respect to tp one obtains the following identities, 



tidyHfa) = idpEi, 

,id v H<t>j) = {Ej - EiX^uidtp^j), i ^ j. 



(19) 



From these one sees that the current matrix elements are 
given by equations, 



la = 2ed v Ei, 

I, j = 2ei(Ei- Ej)Aij, 



or 



^ Sij + i[E,A]ij 



(20) 



(21) 



Thus we conclude that the matrix elements Aij are re- 
lated to the off-diagonal matrix elements of the Josephson 
current operator. 
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Towards the end of this section we present a many 
body Hamiltonian of the junction in the instantaneous 
eigen basis. To this end we define the conjugate momen- 
tum n corresponding to tp, 



dtp 4e 2 



(22) 



and perform a Legendre transformation of the La- 
grangian in Eq. ([£]), then we promote the variables, a i7 etj, 
and (p, n to operators by imposing standard (anti-) com- 
mutation relations to get, 



Uq - 2C 



(2e) 2 I a t s«U 

ij 



(23) 



which the dependence on the time derivative of the phase 
is eliminated from the Hamiltonian. This is achieved 
by using a unitary matrix U(f) satisfying the equation 
id v U = —AU. Computing the derivative and rotating 
back to the original basis we find, 



JSplnHG*- 1 ] 
5x 



= ~^)/>^>> (30) 



where IjUp) is the Josephson current operator defined 
in Eq. (21). Then introducing external current, I e — 



-2ed v U e , we write equation of motion on the form, 



~ip + Tr{Ijf)=I e , Ij 



2e(d v E + i[E,A]) . (31) 



Eqs. ( 28 ) and ( 31 1 together constitute a central technical 



result of this paper. 



III. EQUATION OF MOTION 

Now we perform integration over the fermionic vari- 
ables using the instantaneous eigen basis, 

Z = I XV e^oM+SpM-^ 1 ) =[t><p e ^eftM_ ( 2 4) 



Defining in a standard manner four Green's function 
components, depending on wether the time arguments 
are defined on the forward (a = +) or backward (a = — ) 
part of the contour, 



G(M') = G ab (t,t'), teC a , t 1 GC b 



(25) 



we write Eq. ( 10 1 on the form, 

a (id t - H(<p a ,<p a )^) G ab (t, t') = S ab 5(t - t'). (26) 

Introducing a single particle density matrix through the 
relation, 



we get from Eq. ( 26 1 the Liouville equation, 



if=[H,f], H = E-<pA. 



(27) 



(28) 



A semiclassical dynamical equation for the supercon- 
ducting phase is given by the least action principle for- 
mulated in terms of the Wigner variables, (p a — ip + ax/2, 
and has the form l2l)l. 



SS eS [ip,x\ 



= 0. 



(29) 



To calculate the functional derivative of the fermionic 
part, we perform a rotation to a single particle basis, in 



IV. ADIABATIC LIMIT 

In general, in order to solve the coupled equations for 
the phase plj) and the density matrix ( 28 ) , one needs to 



calculate the static quasiparticle energy spectrum, and 
matrix elements of the interlevel transitions. This is a 
rather difficult task since the latter quantities are com- 
plicated functions of the phase. However, if the quasi- 
particle spectrum has a gap, and the frequency of the 
plasma oscillation is small compared to this gap, in other 
words, if the quasiparticle dynamics is fast on the time 
scale of the phase variation, one can apply an adiabatic 
approximation to find the solution. 

A formal condition for the adiabatic expansion is 
tpAij <C JSj — E j. I n the main approximation, the Hamil- 
tonian in Eq. (28) reads, H = E, and the initial equi- 
librium density matrix, /(0) = f°(E(Q)) defines the so- 
lution that remains constant during the phase e volu tion, 
/(*) = f°- This implies that the trace in Eq. (Ell will 



only contain the diagonal part of the current operator, 
and the Josephson current reduces to the adiabatic form, 



lf(v) - 2eTr(d v £/°) = 2ed v Uj, 
Uj&)=Tr(E(<p)f ). 



(32) 



This equation provides a generalization of the tunnel 
junction equation to the junctions with non-sinusoidal 
current-phase dependence. 

To find the first non-adiabatic correction, it is con- 
venient to expand electronic part in effective action, 



Eq. (24) 



Spln(-iG" 1 ) = Spln(-i(G ad )- 1 ) 

-£-S P (-G a V4) n , (33) 
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where (G ad ) _1 = 6 tj (id t - E t (ip)). The first, adiabatic 
term is given by equation, 

S P ln(-i[G ad ]- 1 ) =-iJdt Uj(ip), (34) 
c 



Let us explicitly evaluate the contribution to Eq. ( 40 ) 



consistent with Eq. (32). To see this, we formally intro- 
duce G| d = G ad (A^?), and rewrite the adiabatic term as 



Spln(-t[G ad ]- 1 ) = J dA^-S P ln(-z(Gf y 1 ) 



(35) 



of the Andreev bound states in a tunnel junction. In 
tunnel junctions, Andreev energy levels are located very 
close to the gap edges [57] having the level spacing, 

e = 2Aw 1 — Z?sin 2 % 2A. The transitions connect 



only Andreev states of the same conducting mode with 
transition matrix elements |38j . 

,Wol „\/T> 

(41) 

Computing the correction to the capacitance using these 
expressions we find the phase dependent correction in the 
zero temperature limit to be, 8C V « (Z?e 2 /32A) cos <p, 
per conducting mode. This is consistent with the result 
of the tunnel model calculation in Refs. l25l l30l. 



Since [G^]u(t,t) = iff docs not change with time by 
virtue of the earlier presented argument, we find, 

Spln(-z[G ad ]- 1 ) 

c ^ ' c 



The first order non-adiabatic term in the series, Eq. (33 1, 
cancels since G ad is diagonal while A is purely off- 
diagonal, which implies that the trace of their product 
is zero. Keeping then only the second order correction 
we find, 



SplnHG- 1 ) = -i j dt Uj(tp) 



c c 

(37) 

When the occupied and unoccupied states are separated 
by a large gap, the product 



G ad (t,i')G a f (t',t) ~ e-' 1 (*«"), 



(38) 

oscillates rapidly on the scale of variations of the phase, 
and we can treat this object in the local approximation. 
This gives us, 



Spln(-iG _1 ) =-ildt 
c 

where 

<5GM = 2e 2 ]T 



\Aij(<p)?Mi-fs) 



(39) 



(40) 



represents a phase dependent correction to the junction 
capacitance. 



V. NON-ADIABATIC DYNAMICS 
A. Linear response 

The non-adiabatic dynamics essentially results from 
the resonant response of low energy quasiparticles to the 
phase variation. In this section we consider the linear 
quasiparticle response and compute the non-adiabatic 
correction to the frequency of Josephson plasma oscil- 
lation. 

Consider small deviations from an equilibrium config- 
uration, ip = ip and f — f° determined by the equation, 
lf(ip ) = Tr(i(<p )p) = I e ((p ). Straightforward lin- 
earization of Eqs. (28 1 and (311 with respect to small 



deviations of the phase, ip(t) — ipo, and the density ma- 
trix, f(t) — f°, leads to the dispersion equation for the 
plasma oscillation, 

(-uj 2 + w 2 + w7o(w)) <p u = 0, (42) 

where 



2e dlf(p ) 
C dip 



(43) 



is the adiabatic plasma frequency, and 7o(w) denotes the 
linear response of the quasiparticles, 



7o M 



4c 2 

~G~ ^ 



Sij — (uj + i0) 



(44) 



The linear response of quasiparticle is a relevant approxi- 
mation at small phase oscillation when the quasiparticles 
have a continuous energy spectrum and the transferred 
energy is dispersed across a large phase space volume re- 
sulting in weak non-equilibrium. As such the dispersion 
equation ( 42 1 can be applied, for example, to the low 



energy itinerant states in the nodal regions of high-Tc 
superconductors, or to broadened Andreev bound states 
in disordered junctions. However, the linear approxima- 
tion does not apply to spectroscopically narrow Andreev 
bound states, whose response is essentially nonlinear even 
at small phase amplitude. 
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B. Resonant interaction with Andreev levels 

Now we consider the nonlinear dynamics of the phase 
driven by small oscillating current I e (t) = I e cos uit, at 
a frequency not far from the resonant frequency S — 
uj — ujp <C 1, in the presence of resonant interaction with 
weakly broadened low energy Andreev levels. Such lev- 
els may exist in transparent electronic conducting modes 
close to tpo = 7T, ha electronic modes with resonant trans- 
missivity, or in surface modes of d-wave superconductors. 
The exact nature of these states does not play any role 
for our analysis. The important properties are: (i) the 
phase variations do not change the electronic momentum 
hence do not induce quasiparticlc transitions among the 
conducting modes, (ii) therefore transitions only occur 
between pairs of Andreev states within the same con- 
ducting mode, (iii) the Andreev levels are well separated 
from the continuum states of the mode. Under these as- 



sumptions, the Hamiltonian in Eq. (28) truncated to the 



Andreev level subspace consists of a sum of independent 
two-level Hamiltonians, and the density matrix factor- 
izes to the product of two-level density matrices param- 
eterized with the conduction mode number, f(n). The 
non-adiabatic current then becomes: 



(It) - I 



de 



2e E(^- f°) -<- A f- + ^ /+) 



(45) 



where, e = E\ — E 2l is the level spacing between two 
Andreev states associated with a specific mode n and, 
iA = A12 , is the corresponding transition matrix element 
(we skip index n for brevity). Similarly f z = /n — f 22 
and /+ = /12 = (/-)*) are the corresponding elements of 
the two-level density matrix satisfying the Bloch-Redfield 
equation, 



/+ = {-ie - T 2 )f+ + 2<pAf z 

f z = -<j>Af% - M*f+ - Ti(/ 2 - /,,o), 



(46) 



where we have added phenomenological decay rates I\ 
and r 2 originating, e.g., from some weak inelastic inter- 
action with the continuum states. 

To separate the fast and slow resonant dynamics, we 
parameterize the phase as, 



¥>(*) = ^r(<Ao(*)e" 



c.c.), 



(47) 



where the complex variable, ^ w (t) = r^e"'' 4 ', depends 
on the amplitude of oscillations, r(t), and the time de- 
pendent phase shift, i)(t). Using a similar separation for 
the fast and slow parts of the off-diagonal elements of the 
density matrix, 



/+(*) = f u (t)e 



(48) 



we get, after expanding to first order in, tp — (po, and 
averaging over fast variables (note Aq — A(tpo) and £0 = 
e(<Po)), 



-i(e - oj - iT 2 )f u - iu>A <Pu>fz 



uj (49) 
fz = i^iAofiuf^ + c.c.) - T 1 (f z - f xfi ). 

The regime relevant for our discussion corresponds to 
slow variation of the phase oscillation envelope, ip u , on 
the time scale of the Andreev state relaxation. Then 
the Andreev state density matrix will adiabatically fol- 
low the evolution of the phase amplitude (in the rotating 
frame), and we restrict ourselves to the quasi-stationary 
solutions, fu)-,fz ~ 0, to find from the first equation in 



(49), 



UlAoPu 



-fz 



(uj + iT 2 ) - £0* 
Inserting this expression into the current we find, 

{ I J )= I f&) + 2eJ2(j^(fz-f°z) 



(50) 



UJE \A \ 2 f z 

(uj + iT 2 ) - £ 



(51) 



Eq. (51 1 illustrates the principal effect of the resonant in- 



teraction between the phase and the Andreev levels: the 
phase oscillation drives the Andreev levels to a nonequi- 
librium state determined by the stationarity condition, 



fz f z 



n 2 (r)(r 2 /r x 



(eo-w^ + ri + fi^rxiyri) 



fl (52) 



Here Q(r) = \u>Aor\ is the amplitude dependent Rabi fre- 
quency of the Andreev two-level system associated with 
specific mode n. The first term inside the bracket in Eq. 



(51) produces a nonlinear modulation of the Josephson 



potential due the nonequilibrium population of the An- 
dreev levels. The second term causes a nonlinear damp- 
ing of the phase oscillation, similar to the imaginary part 
of the linear response, although it now depends on the 
nonequilibrium population of the Andreev levels. 

For the levels close to the resonance, £ « uj, the diag- 
onal elements are approximately given by 



r 2 + n 2 (r) 



,0 

J z 1 



(53) 



where, T = y/TiT 2 . Thus in the limit of O(r) < T, i.e. 
r <C r/|w.4o|, we recover the linear response regime. In 
the opposite limit, fi(r) 3> T, i.e. r ^> T/\ujAo\, the 
levels become saturated, f z ps 0, and can no longer ab- 
sorb energy from the phase oscillation, thus the damping 
decreases for large amplitude of phase oscillation. 
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C. Nonlinear phase dynamics 

To see how the nonlinear quasiparticle response man- 
ifests itself in the junction dynamics we write down 
the equation of motion for the slowly varying ampli- 
tudes, ifu, and introduce a nonlinear response function, 
7(r) = y(r) + ij"(r), defined through the relation, 



2e 
~C 



((Ij)-I?(<p))=ury(u>,r)<p b 



e~ lut + c.c. 



(54) 



in terms of which the averaged equation for the envelope 
becomes, 



2iuj p ip u + [—2u p 6 + io p ^(uj p , r)] ip u 



C 



-J, 



(55) 



The stationary solutions to this equation, <fi u = 0, con- 
nect resonant amplitude and detuning 5, 



6 = h' {r) ± Yr"J {eh/c ^ )2 - {l " [r))2r2 - 



(56) 



The two solutions correspond to the stable/unstable 
branches of the function r(S) as illustrated on Fig. [l] The 
maximum response, r m , is found where the two branches 
coincide, i.e. r m j"(r m ) = eI e /Cuj p . 

To make a quantitative analysis we write = 
J deu(e), where v{e) = J2 n — £ o( n ))- If the density 
of states v{e) is a smooth function close to the resonance 
the integration can be explicitly performed, giving, 



7 (r) = 7o - 



7 "(r) 



r 2 



(57) 



v/r>(r) 2 + r 2 



where j'q is the imaginary part of the linear response 
PI, 



^' = ^4^lM")/,V/2)> 

and bars indicate the values of the functions at the res- 
onance. With this expression we find the maximum re- 
sponse amplitude, 



-1/2 



where 



I, = 



eI R 



Cu p tS 2 7 £ I C 
is the dimensionless driving current, and 

r 



uj p \A \ 



(58) 



(59) 



(60) 



This result shows that the response has an explosive in- 
stability manifested by a divergency of the oscillation am- 
plitude when the driving current amplitude reaches the 



critical value I e — I* . We emphasize that this current 
is much smaller than the Josephson critical current, Jc, 
which sets the scale for the nonlinear behavior of the 
adiabatic junctions. This instability is easy to under- 
stand noticing that the damping produced by the An- 
dreev states decreases with amplitude of oscillation, and, 
on the other hand, it is the damping value that limits 
the resonance response amplitude. To eliminate the di- 
vergency, one has to take into account other damping 
mechanisms, which are weaker than the linear damping 
by the Andreev states. 

If we turn off the external drive, — 0, we find from 



Eq. (55 1 the equation for the decay of the oscillation am- 
plitude, t = -Y'(r)r/2. For r > T/\Ao\u> p , we find that 



the plasma oscillation decays linearly with time with the 
rate, r m — (r^Q / 1 ^4q | ojp ) = const, until it enters the lin- 
ear regime, r < r/|.4o| w p, where the decay crosses over 
to an exponential time dependence, r ~ exp(— y^t). 




FIG. 1: Effect of resonant interaction with spectroscopically 
sharp Andreev bound states on non-linear response of the 
junction. Phase oscillation amplitude as a function of detun- 
ing shown for different amplitudes of driving current. 



VI. LANGEVIN EQUATION 



The classical equation of motion (31 1 is determinis- 



tic and thus does not include the fluctuations originat- 
ing from the coupling of the phase to the quasiparticles. 
In this section we shall outline how these effects can be 
taken into account. The method we adopt results in a 
quantum Langevin equation [26[ [28] |39j , although as we 
shall show, the stochastic force in this case generally has 
non-gaussian properties. 

Expectation values of any function of phase, F(ip), is 
given by 

(F)(h) = J dtp 1 F{<Pi)f>ied(tpi,'Pi,ti), (61) 

where Pred(Vii Vi> ^l) i s the reduced density matrix. 
Noticing that the partition function ^ is given by 



the trace over the reduced density matrix, Z = 
J dipip-cad^i, <pi, t), we are able to write the diagonal 
elements in terms of the Wigner variables ip, x, on the 
form, 



PrediVutfuh) = J V(P J Vxe 
ip(ti)=<pi x(*i)=0 



iS[<p,x] 



(62) 



where the limits on the functional integrals indicate that 
the endpoints, f{t\), x{pi)i °f the trajectories are to be 
held fixed. 

To zeroth order in the saddle point approximation, 
S[ip,x] ~ J dtx(t) (SS[<p, 0]/Sx(t)), only the classical 
path, ip c (t), is realized and the density matrix is writ- 
ten: 



PiBd(<PX,<fl,t) 



J V<p J V X exp (i j dtx{t) 



V>(ti)=V>i 



Vip S 



X(ti)=0 

8S[tp,0] 



sx(t) 



v(ti)=vi 



Sx(t) 



(63) 



where, S[. . .], denotes a delta functional. Average quan- 
tities are then entirely determined by the classical path 
(F)(ti) = F(<p D (ti)). 

To go beyond this deterministic description and in- 
clude fluctuations we can expand the action around the 
saddle point, x = 0, to second order, 



S[<p,x] « / dt X {t) 



5x{t) 



dtdt'xit) 



(64) 



S X (t)5xit 



Here the kernel, 
. 8 2 S[p,0] 

' J2 Tv [G ab (t,t')i{t')G ba (t',t)i{t) 



i2ef 



(65) 



ab 



l 



;Sil<p](t,t'), 



(2e) 

is given by the symmetrized current-current correlation 
function, Si[ip]it,t / ), which is a functional of <p due to 
the dependence of G ab and I on fit). 

We decouple the quadratic term in % by introducing 
an auxiliary variable It which shall later be interpreted 
as a stochastic current [Ml, 



exp 



If 1 

- / dtdt'xit)-, — 
2 J AV ; (2e) 2 



-i/ d «e(*)x(t) 



(66) 



where P[ip,I^] denotes the functional distribution 

P[<p,It] =Af[ (p ] e -hJ^'h(t)S7 1 l'Pm'W) j ( 67 ) 

where N\<p\ = i&etSj 1 ^])" 1 / 2 . The density matrix can 
then be written as 



Predial, fl,t\) 



6S[<p,0] 
~^xW 



hit) 



(68) 



The delta functional selects a single trajectory, tp^, for 
each realization of, 7^, determined by the classical equa- 
tion, 



6S[<p,0] 

Sx(t) 



c_ 

2~e. 



(69) 



where, for the sake of convenience, we assumed an un- 
biased junction. Eq. (69 1, is a stochastic equation and 
averages are given by, 



(70) 



where (. . .)j = fVI^(...)P[(p^,I^]. In contrast to 
the conventional theory of quantum Langevin equations 
the functional distribution, P[ip^,I^\, is in general non- 
gaussian due to the dependence of the symmetrized cur- 
rent correlation function on tp^ = <p[I^\- This is a conse- 
quence of non-equilibrium nature of the fermionic bath 
strongly coupled to the phase variable. 

The stochastic force becomes gaussian under the linear 
response approximation. We consider small deviations 
from a classical equilibrium configuration, (pit) = (po and 
/(*) — /°) an d get the equation, 

(-w 2 +uj 2 p +W70M) Wu) = ( 71 ) 

The functional distribution can be taken at the equilib- 
rium value, P[</?0i£L which then becomes Gaussian and 
the stochastic current, 1^, satisfies the typical relations 
for Gaussian noise: 

(kit)k=0, (I 6 it)hit'))t=Si(t-t'), (72) 
where 

SjV) = coth (^) E i^W - f> 6 ^ - ^) 

= (2e) V coth (i£) \M\fi ~ f!)*S(u, - e l3 ) 

ij 

= Cu coth (£) Im7o (w) . 



(73) 



Thus the fluctuating current is related to the dissipative 
response by the quantum fluctuation dissipation theo- 
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VII. CONCLUSIONS 

We have presented a general theory framework for de- 
scribing non-adiabatic dynamics of Josephson junctions 
with low energy quasiparticle states. The theory applies 
to a wide class of Josephson junctions including trans- 
parent mesoscopic contacts based on 2DEG, nanowires, 
quantum dots, and also junctions of unconventional 
superconductors. It was shown that in the classical 
limit the equation of motion for the phase must be 
solved together with a Liouville equation for density 
matrix of low energy fermionic states. Furthermore, we 



illustrated how the dynamics of such systems can differ 
significantly from the adiabatic (tunnel) junctions, by 
investigating the resonant dynamics of the phase and 
low energy Andreev bound states. It was shown that 
nonlinear, two-state dynamics of the Andreev bound 
states, rather than an adiabatic Josephson energy, de- 
fines the nonlinear macroscopic dynamics of the junction. 

This work was supported by the Swedish Research 
Council (VR), and the European FP7-ICT Project MI- 
DAS. 
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